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We present an approach to calculation of point defect optical and thermal ionization energies based 
on the highly accurate quantum Monte Carlo methods. The use of an inherently many-body 
theory that directly treats electron correlation offers many improvements over the typically-employed 
density functional theory Kohn-Sham description. In particular, the use of quantum Monte Carlo 
methods can help overcome the band gap problem and obviate the need for ad-hoc corrections. We 
demonstrate our approach to the calculation of the optical and thermal ionization energies of the 
F-center defect in magnesium oxide, and obtain excellent agreement with experimental and/or other 
high-accuracy computational results. 



From electronics to optoelectronics to photovoltaics, 
point defects influence and even dominate the proper- 
ties of semiconducting materials jTHB]. Quantitative de- 
scriptions of the effect of point defects on electronic, 
optical, and transport properties is critical to enabling 
point-defect engineering for materials design. How- 
ever, accurate prediction of point-defect energetics, ther- 
mal ionization energies, and optical transition energies 
from first principles remains a challenge. Currently, the 
most widely-used approach based on conventional density 
functional theory (DFT) suffers from poor descriptions of 
band gaps that render difficult the accurate description of 
mid-gap defect states [5[ [7HH]. Here we demonstrate that, 
by contrast, an inherently many-body approach based on 
quantum Monte Carlo (QMC) methods [TUlfTT] can elim- 
inate these problems and enable high-accuracy calcula- 
tions of point defect optical and thermal ionization ener- 
gies. Our computed optical transition energies are in ex- 
cellent agreement with experimental and/or other high- 
accuracy computational results for the same system [T2J , 
and demonstrate that QMC can obtain quantitatively 
accurate descriptions. 

QMC methods comprise a suite of stochastic tools that 
enable calculations of material properties based on the 
many-particle Schrddinger equation. Because of their 
direct treatment of electron correlation, QMC methods 
are among the most accurate electronic structure ap- 
proaches available today, and demonstrate a long and 
distinguished record of ground-breaking and benchmark- 
ing calculations [TUl HTJ |T3] . In comparison to the other 
"beyond-DFT" techniques that are currently explored for 
calculation of point defect properties (DFT+U, hybrid 
DFT, and the GW method), QMC is directly based on 
the true many-body Schrddinger equation and offers the 
possibility of parameter-free accurate band gaps and to- 
tal energies. The application of QMC techniques to point 
defects in solids is still a relatively new field. To date, a 
handful of studies have been carried out to compute de- 
fect formation energies: interstitials in silicon [T4lll6j . va- 
cancies in diamond [TS] , the Schottky defect in MgO [T7] , 



and vacancies and interstitials in aluminum |18j . 

In this letter, we illustrate the application of the QMC 
method to the F-center defect in magnesium oxide (MgO) 
by computing defect formation energies, thermal ioniza- 
tion levels, and optical ionization energies (although the 
approach is general and can be extended to other mate- 
rials of interest). The F-center defect (oxygen vacancy) 
in MgO is a typical example of an intrinsic point de- 
fect in a binary ionic compound [22l427j . Despite its ap- 
parent simplicity, its properties as deduced from optical 
absorption and luminescence studies have proven some- 
what ambiguous. Experimental characterization of the 
F-center in its neutral (F°) and singly ionized (F +1 ) state 
has been complicated by their nearly identical optical ab- 
sorption energies [25H27] . These energies have been cor- 
roborated by recent GW calculations [12]; however, these 
calculations also predicted optical emission energies that 




FIG. 1. (Color online). Left: The electronic band struc- 
ture of a 64 atom MgO supercell containing a single oxygen 
vacancy, calculated within DFT-PBE. The neutral oxygen 
vacancy introduces a localized mid gap defect level of sym- 
metry ai 9 . There is also a triply degenerate excited defect 
level in the conduction band of t\ u symmetry. Right: The 
corresponding ai s and t\ u Kohn-Sham states plotted at the 
T-point, showing the localized nature in the vicinity of the 
vacancy. 
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DFT-PBE 


DMC 


Exp [T9H2T1 


Lattice const (A) 


4.25 


4.22 


4.216 


Coh. En. (cV/MgO) 


9.50 


10.18(5) 


10.5 


OP (eV) 


4.83 


7.96(6) 


7.78 


QP (eV) 




7.9(1) 


7.84 


IP (cV) 




3.28(7) 




EA (eV) 




11.17(7) 





TABLE I. Comparison of lattice constant, atomization en- 
ergy, and band gap in MgO solid according to DFT, DMC, 
and in experiment. All computed results are obtained using 
Ne-core pseudo potentials for magnesium. The optical band 
gap is determined in DFT from the Kohn-Sham levels, and in 
DMC from the extrapolated optical excitation energy. Error 
bars are shown in parenthesis. Note that the IP and EA are 
given with respect to the average potential in the supercell. 



are substantially different from the assigned experimental 
values, causing the authors to suggest a reinvestigation 
of the experimental observations. A particularly com- 
pelling possibility is to explore the F-center defect using 
distinct high-accuracy first-principles techniques to com- 
pare results. Our results - calculated independently using 
QMC methods - corroborate the GW results and further 
invite reassessment of the experimental data for the op- 
tical emission. 

We first compute the properties of the F-center de- 
fect in MgO within a DFT [28j [29] framework as imple- 
mented in the SIESTA code [3D] , employing the Perdew- 
Burke-Ernzerhof |31j approximation to the exchange cor- 
relation potential. The inner core electrons are rep- 
resented by Troullier-Martins pseudopotentials (leaving 
the Mg 3s and O 2s, 2p electrons in valence), and the 
Kohn-Sham orbitals are represented by a linear combi- 
nation of numerical pseudo atomic orbitals expanded in 
a triple-C with polarization Gaussian basis set. For bulk 
rocksalt MgO, in agreement with previous DFT calcu- 
lations [32l [22] we find a lattice parameter of 4.25 A 
(4.22 A in experiment pj5]), an atomization energy of 
9.50 eV/MgO (10.50 eV/MgO in experiment [19 ), and 
a direct band gap of 4.83 eV at the T-point (a consider- 
able underestimate of the experimental band gap of 7.78 
eV [H]). 

The DFT band gap underestimate has a severe con- 
sequence on the prediction of mid-gap defect states, de- 
fect energetics (particularly for occupied defect levels), 
and defect-induced optical absorption and emission en- 
ergies. Broadly, the band gap underestimate arises be- 
cause in typically-employed mappings of the interact- 
ing many-body Schrodinger equation to the DFT single- 
particle "effective" -potential Kohn-Sham equations, each 
electron also interacts with itself (self-interaction er- 
ror [34l438| ). This results in an extraneous Coulomb 
repulsion that overly delocalizes electronic states. The 
self-interaction error, in addition to the absence of a 
derivative discontinuity in the exchange-correlation po- 
tential [51H55] , results in underestimated band gaps that 



have plagued DFT calculations. 

In Fig. [I] we show the DFT-computed electronic band 
structure of a 63 atom MgO supercell containing an F°- 
center defect (neutral oxygen vacancy). In agreement 
with previous DFT calculations [12], the F°-center intro- 
duces a fully-occupied mid-gap defect level of a\ g sym- 
metry into the electronic band structure; higher in the 
conduction band we also find a triply-degenerate excited 
defect level of t\ u symmetry. 

The QMC calculations reported here are computed 
within fixed node diffusion Monte Carlo (DMC) as imple- 
mented in the QWalk code [35], with single-determinant 
Slater-Jastrow trial wave functions constructed from the 
DFT orbitals, variance-minimized Jastrow coefficients, 
and a time step of 0.01 au. To establish that our choice 
of pseudopotentials is reasonable, we first calculated the 
bond length, electron affinity, and binding energy of the 
MgO molecule within DMC. We tested both Ne and He 
-core pseudopotentials for the Mg atom in the molecule, 
and found that (although both give good results) the 
small core pseudopotential gives a somewhat better de- 
scription (see supplementary information). This suggests 
that including the Mg 2s and 2p electrons improves the 
description slightly; however, for the solid, the compu- 
tational cost of the He-core pseudopotential was pro- 
hibitive. 

To test the properties of the non-defective MgO solid 
in QMC, we calculate the atomization energy, the opti- 
cal band gap, and the quasiparticle band gap. All en- 
ergies are calculated using the extrapolation framework 
described in Refs dUHH] with supercells containing 16, 
32, and 64 atoms. Further details are provided in the 
supplementary information. We calculate the following 
states: the ground state E g , the T-point optically ex- 
cited state Ep^r, the positively charged state E+, and 
the negatively charged state E_. The ionization poten- 
tial (IP), electron affinity (EA), quasiparticle gap (QP), 
and optical gap (OP) are given by 

IP = E g - E + (1) 

EA = E- - E g (2) 

QP = EA - IP (3) 

OP = E r ^r - E g (4) 

The results are summarized in Table [I] and show ex- 
cellent agreement overall with the experimental values. 
The slight underestimation of the atomization energy is 
likely due to the Ne-core pseudopotential, since the MgO 
molecule showed a similar effect (see supplementary in- 
formation) , while the gap calculations are close to exper- 
iment, overestimating the gap slightly. 

Now we turn to the F-center defect. We again use 
64 (perfect) and 63 (with F-center) atom supercells, and 
compute the defect formation energy AEjj g according to 

AE D<q = (E D<q - E perf ) - nm + q(E v + E F ) (5) 
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Mg-Mg (A) 


O-O (A) 


Relaxation (eV) 


perfect 


2.98 


5.96 




F 


2.99 


5.96 


0.003 


F+ 


3.09 


5.90 


0.545 




3.17 


5.84 


1.182 



TABLE II. DFT-computed lattice relaxations for the F, F+, 
and F +2 center. The Mg-Mg distance denotes the separa- 
tion between Mg atoms that neighbor the missing O atom; 
similarly the O-O distance denotes the separation between O 
atoms that neighbor the missing O atom. 




FIG. 2. (Color online). Comparison of F-center defect for- 
mation energies and thermal ionization energies in MgO com- 
puted in DFT and DMC. The domain of the Fermi energy (x- 
axis) is determined by the band gap of the system according 
to the computational framework; clearly the DMC descrip- 
tion of the gap is better and obviates the need for "band-gap 
corrections". In comparison to DFT, DMC modifies some- 
what the absolute value of the defect formation energies, but 
maintains thermal ionization levels near mid gap. (Note that 
the error bars of the DMC-computed formation energies are 
smaller than the line widths in (b).) 



where Er>, q is the (computed) total energy of the super- 
cell containing a defect D in the charge state q, E per f is 
the (computed) total energy of the perfect supercell, and 
rii is the number of atoms of species i added to (rii > 0) 
or removed from (m < 0) the supercell to create the de- 
fect Different environmental conditions are accom- 
modated by the set of chemical potentials /ij for each 
element by assuming that each is in equilibrium with a 
physical reservoir such as a gas or a bulk phase. Ey is 
the energy of the valence band maximum (the ionization 
potential in DMC), and Ep is the Fermi energy refer- 



enced to Ey so that < Ep < E g where E g is the band 
gap. 

The thermal ionization energies, which determine the 
shallow or deep nature of a defect, correspond to the 
Fermi energies at which the energetically most favored 
charge state of the defect changes. According to our 
DFT calculations, the creation of an F°-center results 
in the formation of a filled mid-gap defect level (shown 
m Fig. [IJ. There is very little lattice relaxation that 
takes place upon removal of the O, as indicated in Ta- 
ble [IT] However, when an electron is removed from the 
supercell to form the F +1 -center, in DFT we find a large 
lattice relaxation as the positively charged Mg ions move 
outwards away from, and the negatively charged O ions 
move inwards towards, the positively charged vacancy 
in conjunction with a 0.55 eV drop in energy. Further 
ionizing the defect into the F +2 state in DFT results in 
further lattice relaxations accompanied by an energy re- 
covery of 1.18 eV. The DFT defect formation energies 
obtained from Eq. [5] are plotted in Fig. [2^,, showing ther- 
mal ionization levels near the middle of the gap. These 
formation energies are computed in the Mg-rich limit so 
that (JiMg is given by the chemical potential of solid ele- 
mental magnesium (and [iMg + Ho = HMgO 

). In Fig.|k 

we have used the as-computed DFT band gap for MgO, 
without correction schemes to artificially open the gap. 
Note that we have also ignored the interaction energy 
between charged defects arising from the use of periodic 
boundary conditions, which would shift upwards the lines 
for the formation energy of the F +1 and F +2 center (and 
reduce the defect transition levels). We will ignore this 
interaction for the QMC calculations as well, and com- 
pare both methods on an even footing [44] . 

For the QMC calculations of total energies of per- 
fect and defected supercells (E per /, ^D,q) and the cor- 
responding thermal ionization levels, we use the relaxed 
lattice geometries and orbitals obtained from the DFT 
calculations. For computational efficiency, only the T 
point Kohn-Sham orbitals are used, rather than twist- 
averaging, since in Eq. [5] the errors in the difference 
{Ed. q~ E per f) largely cancel. For instance, in DFT these 
errors are 0.26 eV for the F°-center, and < 0.1 eV for 
the F +1 and F +2 -center, and do not substantively alter 
our results. Figure [2}d shows the defect formation ener- 
gies and the thermal ionization levels as computed within 
QMC. The largest single difference between the DFT and 
the QMC results is, of course, the domain for the Fermi 
energy < Ep < E g . For this material, both DFT and 
QMC place the transition levels near mid-gap |45j . 

From Fig. [2] the overall formation energy of the neu- 
tral defect (F°) is higher in QMC by approximately 1.5 
eV. This increase in formation energy may be largely 
attributed to the fact that DFT-PBE renders the MgO 
system more delocalized while QMC captures the true 
ionic nature of the solid. The overly-delocalized (more 
metal-like) description in DFT makes the penalty for 
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FIG. 3. (Color online). Optical absorption and emission 
energies (in eV) computed in DMC for the F-center defect in 
MgO. The DMC absorption energies are in excellent agree- 
ment with experiment and recently published GW results. 
The DMC emission energies are in disagreement with the ex- 
perimentally assigned values, but match closely to the GW 
results. 

bond breaking too small, and consequently the defect for- 
mation energy too small as well. This over-delocalization 
results in the underestimated DFT band gap; the effect 
on computed energetics are significant for defects with 
occupied mid-gap defect levels such as the F° center. As 
a result of the underestimated band gap, in DFT the 
Kohn-Sham level of the deep, doubly occupied F° level is 
squeezed too close to the valence band maximum, directly 
resulting in a low calculated formation energy. Our find- 
ing here is similar to the findings for the formation ener- 
gies of neutral, interstitial Si atoms in silicon [T4HT6] , for 
which DMC calculations show that DFT underestimates 
formation energies in the case of occupied mid-gap defect 
levels. This analysis is also consistent with recent QMC 
results for defect formation energies in aluminum (TS] - 
for metallic systems, for which DFT derealization prob- 
lems are less significant, the DMC results are more closely 
matched. 

For the charged defects F +1 and F +2 , the compari- 
son between the DFT and QMC results is more complex. 
The difference between the QMC and DFT -computed 
formation energies arises from an interplay between the 
difference in predicted band gap, the occupation of the 
Kohn-Sham defect level (zero, one, and two electrons for 
the F +2 , F +1 , and F° center respectively), and the fact 
that the charged-defect interactions arising from the use 



of periodic boundary conditions are likely different for 
the two schemes. That is, since the degree of localiza- 
tion/delocalization is different, screening is also expected 
to be different for the two methodologies - most likely en- 
hanced in the case of DFT. 

We now turn to the QMC description of the opti- 
cal ionization energies, corresponding to vertical Franck- 
Condon transitions on a configuration coordinate dia- 
gram as illustrated in Fig. [3j An optical transition oc- 
curs when a photon is absorbed or emitted by the defect; 
because this transition essentially takes place instanta- 
neously on the scale of lattice relaxations it occurs at 
fixed atomic coordinates (and hence is represented as a 
vertical transition). Such a transition places the system 
in an excited vibrational state; for example, F°-center ab- 
sorption illustrated in Fig. [3] refers to the absorption of a 
photon and the promotion of an electron from the filled 
mid-gap level to the conduction band, leaving behind an 
electron in the conduction band and an F +1 center in an 
excited vibrational state (which soon decays to the F +1 
vibrational ground state). Therefore, we compute the op- 
tical transitions by using the relaxed coordinates of the 
initial state, and occupying the Kohn-Sham orbitals as 
appropriate to describe the vibronic state. 

Our DMC absorption energies (Fig[3| are in excellent 
agreement with experiment and remarkably close to the 
GW-computed values, demonstrating the high-accuracy 
potential of the DMC methodology. The DMC emission 
energies are also remarkably similar to the GW-computed 
values, but in disagreement with the experimental num- 
bers. The disagreement between the GW and experi- 
mental values led the authors in [5] to suggest that the 
low energy signal around 2.3-2.4 eV that is observed in 
fact arises when electrons in the defect level recombine 
with holes in the valence band. We find it notable that 
two distinct many-body approaches (namely QMC and 
GW) have yielded similar results for the optical emission 
transitions in question. 

This leads us to suggest two possibilities. First, we find 
it likely that, as suggested by the authors in Ref. jS], the 
original emission-peak assignment should be revisited. A 
second possibility is based on the fact that both our QMC 
and the GW results are built from DFT-relaxed atomic 
geometries (Table [n]). It is possible that the GW and 
QMC results compare favorably because both methods 
are using similar DFT-relaxed lattice geometries. If the 
relaxations are not properly described in DFT, then the 
many-body energies may be similar but incorrect. How- 
ever, the possibility that the lattice geometries are prob- 
lematic seems unlikely given the exceptional agreement 
with experiment for the absorption transitions. 

In conclusion, we demonstrate the application of quan- 
tum Monte Carlo methods to the calculation of the ther- 
mal and optical ionization energies of point defects in 
solids. The striking agreement between two highly accu- 
rate methods, quantum Monte Carlo and GW, suggests 



that predictive calculations of point defect properties are 
now in reach. Due to its inherently many-body approach 
and accurate treatment of electron correlation, quantum 
Monte Carlo shows large promise for the quantitative 
first-principles calculation of point defect properties. 
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Supplementary Information 

I. Properties of MgO Molecule as Described in 
Diffusion Monte Carlo 

We compute the properties of the MgO molecule in 
diffusion Monte Carlo (DMC) using both small (he- 
lium) and large (neon) -core pseudopsotentials; the re- 
sults are given in Table |III| While both pseudopoten- 
tials give excellent agreement to the experimental bond 
length and electron affinity, we find some difference for 
the binding energy and the ionization energy. Using a 
He-core, rather than Ne-core, pseudopotential for Mg in- 
creases the molecular binding energy from 2.28±0.01 eV 
to 2.43±0.01 eV, in comparison to the experimental value 
of 2.54±0.22 eV [HIMgS]. This suggests that allowing 
the Mg 2s and 2p electrons to participate in the bond- 
ing allows more recovery of the binding energy. We also 
find that using the Ne-core pseudopotential introduces a 
small 0.04(2) eV error in the ionization energy compared 
to the He core pseudopotential. 





bond 
length (A) 


binding 
energy (eV) 


electron 
affinity (eV) 


DMC, Ne-core PP 


1.75 


2.28(1) 


1.76(1)) 


DMC, Ar-core PP 


1.75 


2.43(1) 


1.72(1) 


Exp. (Refs. HHSHS]) 


1.75 


2.56(21) 


1.630(25) 



TABLE III. Comparison of bond length, binding energy, and 
electron affinity for the MgO molecule according to DMC and 
in experiment. Two sets of DMC results are provided, corre- 
sponding to the use of neon (large) and helium (small) -core 
pseudopotentials. Error bars are shown in parenthesis. 



II. Properties of Crystalline MgO as Described 
in Diffusion Monte Carlo 

The extrapolation framework described in Refs. |40l - 
|4"2"] is used to compute the atomization energy for the 
MgO solid, as shown in Fig. |4j We use supercells con- 
taining 16, 32, and 64 atoms using twist-averaged bound- 
ary conditions. The dependence of the binding energy 
on the supercell size reflects the spurious electron cor- 
relation that appears in many-body theories when pe- 
riodic boundary conditions are applied. This spurious 
correlation disappears in the infinite size supercell limit. 
The results presented in Fig. [4] are computed using twist- 
averaged boundary conditions. The extrapolated value of 
the atomization energy in DMC is 10.18±0.05 eV per for- 
mula unit, in comparison to the experimental and DFT 
values of 10.5 eV and 9.48 eV, respectively [TM2T], Al- 
though DMC improves the atomization energy in com- 
parison to DFT, it is most likely necessary to include the 
Mg 2s and 2p electrons in valence to obtain atomization 
energies closer to experiment. 

The optical gap, ionization potential, and electron 
affinity are similarly computed by extrapolation, as illus- 



trated in Fig. [5] Here we use only the T-point orbitals to 
construct the many-particle wave function, to save com- 
putational cost and since the MgO solid exhibits a direct 
band gap at the T-point. For these quantities, additional 
finite size effects are present including (1) periodic im- 
age interactions between the electron-hole pair for the 
optical gap and (2) the electrostatic interaction between 
charged supercells in the calculation of the ionization po- 
tential and electron affinity. Extrapolating to the infinite 
supercell limit, we obtain an optical gap of 7.96±0.06 eV 
and a quasiparticle gap of 7.89±0.10 eV, in close agree- 
ment with the experimental band gap of 7.8 eV. 




1/64 1/32 1/16 
1/(number of atoms) 



FIG. 4. (Color online). The extrapolated value of the at- 
omization energy as computed in DMC is 10.18±0.05 eV, in 
comparison to the experimental and DFT values of 10.5 eV 
and 9.48 eV, respectively [19H2T] . 



> 

>. 
CO 

03 
LU 



b) 



8.0 
7.5 
7.0 
6.5 
6.0 



r-point optical excitation 





12 




10 


> 










8 




CD 




CD 


6 


C 




LU 






4 




2 



1/128 1/64 1/32 1/16 
r-point quasiparticle gap 







E(N+1)-E(N) 






E(N)-E(N-1) 






1/64 


1/32 


1/16 



1 /(number of atoms) 

FIG. 5. (Color online). Optical and quasiparticle gap 
in MgO obtained by extrapolating supercells to infinite 
size. The data points for finite-size supercells are com- 
puted with DMC. The extrapolated value of the optical 
gap is 7.96±0.06 eV. The extrapolated ionization poten- 
tial IP=E(N)-E(N-1)=3.28±0.07 eV. The extrapolated elec- 
tron affinity EA=E(N+1)-E(N)=11.17±0.07 eV. This gives a 
quasiparticle gap of QP=EA-IP=7.89±0.f eV. (Note the EA 
and IP are referenced to the average electrostatic potential in 
the supercell.) 
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